! (C) Copyright 2005- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

MODULE GET_POINTS_MOD

CONTAINS

SUBROUTINE GET_POINTS(INFO_FILE,NAM_FILE,NP_TOT,NPOINTS,POINTS,& 
                    & NP_LAND_OUT,NP_LAKE_OUT,NP_LAKE_FRAC_OUT,NP_OCEAN_OUT)
USE NETCDF
IMPLICIT NONE

CHARACTER(LEN=*),INTENT(IN) :: INFO_FILE,NAM_FILE
INTEGER,INTENT(OUT) ::NP_TOT,NPOINTS
INTEGER,DIMENSION(:),ALLOCATABLE,INTENT(OUT) :: POINTS
INTEGER,INTENT(OUT),OPTIONAL :: NP_LAND_OUT,NP_LAKE_OUT,NP_LAKE_FRAC_OUT,NP_OCEAN_OUT

INTEGER :: NCID,VARID,DIMID
REAL(KIND=8),DIMENSION(:),ALLOCATABLE :: LSM,CLAKE,LAT_TOT,LON_TOT
INTEGER :: NP_LAND,NP_LAKE,NP_OCEAN,NP_LAKE_FRAC,JP,JPP
LOGICAL,DIMENSION(:),ALLOCATABLE :: LDLAND,LDLAKE,LDOCEAN,LDPOINT

LOGICAL      :: LLAND,LLAKE,LOCEAN
REAL(KIND=8) :: RFRACMIN_LAKE,RLATMIN,RLATMAX,RLONMIN,RLONMAX

NAMELIST /NAMGP/LLAND,LLAKE,LOCEAN,RFRACMIN_LAKE,RLATMIN,RLATMAX,RLONMIN,RLONMAX        


!* SET DEFAULT VALUES FOR NAMGP 
LLAND=.TRUE.
LLAKE=.TRUE.
LOCEAN=.TRUE.
RFRACMIN_LAKE=0.
RLATMIN=-100.
RLATMAX=100.
RLONMIN=-1.
RLONMAX=370. 

OPEN(99,FILE=TRIM(NAM_FILE),STATUS='OLD')
READ(99,NAMGP)
CLOSE(99)

IF ( RLONMAX< 0. .OR. RLONMIN < 0. .AND. RLONMIN .NE. -1 ) THEN
  PRINT*,'WRONG LONGITUDE INTERVAL !'
  PRINT*,'ONLY BETWEEN 0/360'
  STOP 9 
ENDIF
IF ( .NOT. LLAKE .AND. RFRACMIN_LAKE > 0. ) THEN
  PRINT*,'LLAKE = FALSE AND RFRACMIN_LAKE > 0. !!!!'
  STOP 9
ENDIF

!***************************************
! OPEN FILE 
CALL NCERROR( NF90_OPEN(INFO_FILE,NF90_NOWRITE,NCID) )

! GET LSM AND CLAKE LAT AND LON TO CALCULATE NUMBER OF POINTS 
CALL NCERROR( NF90_INQ_DIMID(NCID,'x',DIMID) )
CALL NCERROR( NF90_INQUIRE_DIMENSION(NCID,DIMID,LEN=NP_TOT) )
ALLOCATE(LSM(NP_TOT))
ALLOCATE(CLAKE(NP_TOT))
ALLOCATE(LAT_TOT(NP_TOT))
ALLOCATE(LON_TOT(NP_TOT))
CALL NCERROR( NF90_INQ_VARID(NCID,'lsm',VARID) )
CALL NCERROR( NF90_GET_VAR(NCID,VARID,LSM) )
CALL NCERROR( NF90_INQ_VARID(NCID,'cl',VARID) )
CALL NCERROR( NF90_GET_VAR(NCID,VARID,CLAKE) )
CALL NCERROR( NF90_INQ_VARID(NCID,'lat',VARID) )
CALL NCERROR( NF90_GET_VAR(NCID,VARID,LAT_TOT) )
CALL NCERROR( NF90_INQ_VARID(NCID,'lon',VARID) )
CALL NCERROR( NF90_GET_VAR(NCID,VARID,LON_TOT) )


ALLOCATE(LDLAND(NP_TOT))
ALLOCATE(LDLAKE(NP_TOT))
ALLOCATE(LDOCEAN(NP_TOT))
ALLOCATE(LDPOINT(NP_TOT))

!* SEARCH FOR THE POINTS IN THE CONDITIONS ! 
NP_LAND=0
NP_LAKE=0
NP_LAKE_FRAC=0
NP_OCEAN=0
NPOINTS=0
DO JP=1,NP_TOT
  LDLAND(JP)=LSM(JP) >= 0.5 
  LDLAKE(JP)=CLAKE(JP) >= 0.5*(1.-LSM(JP))
  LDOCEAN(JP)=LSM(JP) < 0.5 .AND. (.NOT. LDLAKE(JP) .OR. .NOT. LLAKE )
  LDPOINT(JP)=.FALSE.
  
  IF ((LAT_TOT(JP) <= RLATMAX .AND. LAT_TOT(JP) >= RLATMIN).AND. &
 &    ((RLONMAX>RLONMIN.AND.LON_TOT(JP)<=RLONMAX.AND.LON_TOT(JP)>=RLONMIN).OR. &
 &     (RLONMAX<RLONMIN.AND.(LON_TOT(JP)>= RLONMAX.OR.LON_TOT(JP)<=RLONMIN)))) THEN
    IF ( CLAKE(JP) >=RFRACMIN_LAKE ) THEN 
      IF ( LDLAND(JP) ) THEN    
        IF ( LLAND ) THEN
          NP_LAND=NP_LAND+1
          LDPOINT(JP)=.TRUE.
      IF ( LDLAKE(JP) ) NP_LAKE_FRAC=NP_LAKE_FRAC+1
        ENDIF
      ELSEIF( LDLAKE(JP) .AND. LLAKE ) THEN
        NP_LAKE=NP_LAKE+1
        LDPOINT(JP)=.TRUE.   
      ELSEIF( LDOCEAN(JP) .AND. LOCEAN) THEN
        NP_OCEAN=NP_OCEAN+1
        LDPOINT(JP)=.TRUE.
      ENDIF
    ENDIF
  ENDIF
  IF ( LDPOINT(JP) ) NPOINTS=NPOINTS+1
ENDDO
  

ALLOCATE(POINTS(NPOINTS))
JPP=0
DO JP=1,NP_TOT
  IF ( LDPOINT(JP) ) THEN
    JPP=JPP+1
    POINTS(JPP)=JP
!     WRITE(*,'(I6,1X,4(F9.4))'),JPP,LAT_TOT(JP),LON_TOT(JP),LSM(JP),CLAKE(JP)
  ENDIF
ENDDO

IF ( PRESENT(NP_LAND_OUT) ) NP_LAND_OUT=NP_LAND
IF ( PRESENT(NP_LAKE_OUT) ) NP_LAKE_OUT=NP_LAKE
IF ( PRESENT(NP_LAKE_FRAC_OUT) ) NP_LAKE_FRAC_OUT=NP_LAKE_FRAC
IF ( PRESENT(NP_OCEAN_OUT) ) NP_OCEAN_OUT=NP_OCEAN


!******************************************
! CLOSE FILE 
CALL NCERROR( NF90_CLOSE(NCID) )





END SUBROUTINE GET_POINTS 




SUBROUTINE NCERROR(STATUS,STRING)
USE NETCDF
IMPLICIT NONE

INTEGER,INTENT(IN) :: STATUS
CHARACTER(LEN=*),INTENT(IN),OPTIONAL :: STRING 

IF ( STATUS /= 0 ) THEN 
  PRINT *, TRIM(NF90_STRERROR(STATUS))
  IF( PRESENT(STRING) ) PRINT*,TRIM(STRING)
  STOP
ENDIF

END SUBROUTINE NCERROR



END MODULE GET_POINTS_MOD
